library(osmdata)
library(opentripplanner)
library(sf)
library(tidyverse)
library(ggthemes)
library(ggspatial)
library(viridis)
library(dplyr)Tstops <- st_read("C:\\Users\\emmac\\Documents\\GitHub\\emmacolley-vis\\T_stops\\MBTA_Rapid_Transit.shp")affordable_2020 <- st_read("C:\\Users\\emmac\\Documents\\GitHub\\emmacolley-vis\\affordable_2020\\massbuilds-shp-20201004-a44e92.shp") %>%
filter(MUNICIPAL == "Boston") %>%
st_transform(crs = "+proj=longlat +datum=WGS84") #project to wgs84
MA_state_plane <- "+proj=lcc +lat_1=41.71666666666667 +lat_2=42.68333333333333 +lat_0=41 +lon_0=-71.5 +x_0=200000 +y_0=750000 +ellps=GRS80 +units=m +no_defs"opq(bbox = 'Boston MA USA') %>%
add_osm_feature(key = 'highway') %>%
osmdata_xml(file = 'OTP/graphs/default/boston_streets.osm')## The OTP will be saved to OTP/otp.jar
path_data <- file.path(getwd(), "OTP")
path_otp <- paste(path_data, "otp.jar",sep = "/")
otp_build_graph(otp = path_otp, dir = path_data, memory = 1024)## 2020-10-07 01:09:58 Basic checks completed, building graph, this may take a few minutes
## The graph will be saved to C:/Users/emmac/Documents/GitHub/emmacolley-vis/OTP
## 2020-10-07 01:10:12 Graph built
## 2020-10-07 01:10:13 OTP is loading and may take a while to be useable
## Router http://localhost:8080/otp/routers/default exists
## 2020-10-07 01:11:14 OTP is ready to use Go to localhost:8080 in your browser to view the OTP
## Router http://localhost:8080/otp/routers/default exists
iso_5min_walk <-
otp_isochrone(otpcon = otpcon, fromPlace = affordable_2020,
mode = "WALK", cutoffSec = 300) %>%
st_transform(crs = MA_state_plane) %>%
mutate(mode = "walk")## Warning in otp_isochrone(otpcon = otpcon, fromPlace =
## affordable_2020, mode = "WALK", : Failed to get isochrone
## with error: {"type":"FeatureCollection","features":
## [{"type":"Feature","geometry":null,"properties":{"time":
## 300},"id":"fid-24cdf481_175016df537_-7fac"}]}
iso_5min_drive <-
otp_isochrone(otpcon = otpcon, fromPlace = affordable_2020,
mode = "CAR", cutoffSec = 300) %>%
st_transform(crs = MA_state_plane) %>%
mutate(mode = "drive")
iso_all_modes <- rbind(iso_5min_drive, iso_5min_walk)right_side <- st_bbox(iso_all_modes)$xmax
left_side <- st_bbox(iso_all_modes)$xmin
top_side <- st_bbox(iso_all_modes)$ymax
bottom_side <- st_bbox(iso_all_modes)$ymin
ggplot(iso_all_modes) +
annotation_map_tile(zoomin = 0, type = "cartolight", progress = "none") +
geom_sf(aes(fill = mode), alpha = 0.5) +
geom_sf(data = affordable_2020) +
geom_sf(data = Tstops, color = "lightslategrey") +
coord_sf(xlim = c(left_side, right_side),
ylim = c(bottom_side, top_side), expand = FALSE) +
scale_fill_manual(values = c("skyblue", "royalblue"), name = "Area that is reachable within", labels = c("5 min by car", "5 min by foot")) +
theme_map() +
labs(caption = "Basemap Copyright OpenStreetMap contributors")## Loading required namespace: raster
I’ve labeled each housing development with their neighborhood to begin to understand differences between neighborhoods and commuting distances.
iso_areas <- iso_all_modes %>%
mutate(area = st_area(iso_all_modes)) %>%
st_set_geometry(NULL) %>%
pivot_wider(names_from = mode, values_from = area)
ggplot(iso_areas,
aes(x = as.numeric(walk), y = as.numeric(drive))) +
geom_point(color = "lightslategrey", size = 3) +
geom_text(aes(label = c("Allston", "South Boston", "Dorchester", "Dorchester", "5", "6", "7","8","9","10")), hjust=-.3) +
scale_x_continuous(name = "Area within a ten-minute walking distance\nof housing development\n(square km)",
breaks = breaks <- seq(10000, 130000, by = 20000),
labels = breaks / 1000000) +
scale_y_continuous(name =
"Area within a five-minute driving distance\nof housing development\n(square km)",
breaks = breaks <- seq(0, 700000, by = 100000),
labels = breaks / 1000000) +
theme_minimal()## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_text).
## Simple feature collection with 19 features and 4 fields
## geometry type: GEOMETRY
## dimension: XY
## bbox: xmin: 228593.8 ymin: 893417.4 xmax: 237829.8 ymax: 901383.4
## CRS: +proj=lcc +lat_1=41.71666666666667 +lat_2=42.68333333333333 +lat_0=41 +lon_0=-71.5 +x_0=200000 +y_0=750000 +ellps=GRS80 +units=m +no_defs
## First 10 features:
## id time fromPlace mode geometry
## 1 1 300 42.3165603,-71.0966731 drive MULTIPOLYGON (((232967.4 89...
## 2 1 300 42.3319784,-71.0502934 drive MULTIPOLYGON (((237064 8977...
## 3 1 300 42.3277973,-71.0563240 drive MULTIPOLYGON (((236345.8 89...
## 4 1 300 42.3049051,-71.0784474 drive MULTIPOLYGON (((234392.1 89...
## 5 1 300 42.3358,-71.1094 drive MULTIPOLYGON (((232188 8987...
## 6 1 300 42.360071,-71.149450 drive MULTIPOLYGON (((228881.4 90...
## 7 1 300 42.2974755,-71.1155956 drive MULTIPOLYGON (((231478 8935...
## 8 1 300 42.3007916,-71.1119202 drive MULTIPOLYGON (((231631.9 89...
## 9 1 300 42.3441666,-71.0636718 drive MULTIPOLYGON (((235949.3 89...
## 10 1 300 42.3387177,-71.0697751 drive MULTIPOLYGON (((235453 8984...
## although coordinates are longitude/latitude, st_covers assumes that they are planar
ggplot(iso_transformed) +
annotation_map_tile(zoomin = 1, type = "cartolight", progress = "none") +
geom_sf(aes(fill=transit_score), alpha=.5) +
theme_map() ## although coordinates are longitude/latitude, st_intersects assumes that they are planar
## although coordinates are longitude/latitude, st_intersects assumes that they are planar
iso_transformed <- iso_transformed %>%
mutate(num_Tstops = lengths(st_covers(iso_transformed, Tstops))) %>%
mutate(has_Tstops = num_Tstops > 0)## although coordinates are longitude/latitude, st_covers assumes that they are planar
## [1] 6
ggplot(iso_transformed) +
annotation_map_tile(zoomin = 1, type = "cartolight", progress = "none") +
geom_sf(data = iso_transformed,
aes(fill = has_Tstops)) +
scale_fill_manual(values = c("gray85", "darkblue"),
name = "Boston Neighborhoods\nby presence of a community center",
labels = c("Without a T stop",
"With a T stop")) +
theme_map()## although coordinates are longitude/latitude, st_intersects assumes that they are planar
## although coordinates are longitude/latitude, st_intersects assumes that they are planar
## [1] "SUCCESS: The process \"java.exe\" with PID 21540 has been terminated."